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Abstract 

Mean-field theory is a powerful tool for studying large neural networks. However, when the system 
is composed of a few neurons, macroscopic differences between the mean-field approximation and the 
real behavior of the network can arise. Here we introduce a study of the dynamics of a small firing-rate 
network with excitatory and inhibitory populations, in terms of local and global bifurcations of the neural 
activity. Our approach is analytically tractable in many respects, and sheds new light on the finite-size 
effects of the system. In particular, we focus on the formation of multiple branching solutions of the 
neural equations through spontaneous symmetry-breaking, since this phenomenon increases considerably 
the complexity of the dynamical behavior of the network. For these reasons, branching points may reveal 
important mechanisms through which neurons interact and process information, which are not accounted 
for by the mean-field approximation. 


Author Summary 

The mesoscopic scale represents the bridge between microscopic neural activity and the highest cognitive 
functions that emerge at the macroscopic level in the brain. For this reason, understanding the dy¬ 
namics of small neural networks at this intermediate scale of organization is of fundamental importance. 
However, counter-intuitively, the behavior of small neural networks can be much more difficult to study 
mathematically than that of large networks. Here we introduce a finite-size firing-rate model and we 
study its local and global bifurcations by a combination of numerical and analytical techniques. This 
analysis shows the formation of strong and previously unexplored finite-size effects, that are particularly 
hard to detect in large networks. Their study advances the state of the art in the comprehension of neural 
circuits, beyond the results provided by the mean-field approximation and the techniques introduced so 
far for the quantification of finite-size effects. 

1 Introduction 

The structural complexity of the brain is reflected by its organization at multiple spatial scales [^. At 
the highest level, the brain can be divided in macroscopic areas containing millions of neurons. These 
areas accomplish very complex roles, ranging from motor control to cognitive functions. On the opposite 
side, namely at the lowest level, the brain is made of its elementary blocks, the neurons, which provide 
a description of this complex organ at the microscopic scale. The level of spatial organization that lies 
between these macroscopic and microscopic pictures of the brain is called mesoscopic scale, and in the last 
years an increasing number of studies (e.g. [2j|^) have recognized its importance for the comprehension 
of the brain. This is the spatial scale that links the macroscopic and microscopic scales, and therefore 
that may explain how the highest cognitive functions arise from the cooperation and the exchange of 
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information between neurons. For this reason, finding an appropriate mathematical description of the 
brain at the mesoscopic scale is of fundamental importance for unveiling its emergent properties. 

At the mesoscopic scale, the brain is often described as a collection of neural masses, i.e. homogeneous 
neuronal populations within a cortical column [^. Usually, these groups of neurons are described by the 
so called neural-mass models [^. A typical example is the well-known Jansen-Rit model [9 -11 , which 
describes a cortical minicolumn by a population of excitatory pyramidal cells, that receive inhibitory 
and excitatory feedback from local interneurons, and excitatory input from neighboring areas such as the 
thalamus. This model was originally intended to describe the formation of alpha rhythms within a cortical 
minicolumn of the visual cortex through the interaction of the excitatory and inhibitory populations, but 
provides only a coarse description of their dynamics. In |12| the authors proposed an improvement of 
the Jansen-Rit model based on the mean-held theory, which may provide a better understanding of 
MEG/EEG cortical recordings. The mean-held theory is a mathematical tool that approximates the 
behavior of large networks 13-15 , and proves very convenient from the computational point of view 


since it describes the activity of large populations by means of a reduced number of equations, which 
are usually simpler to solve than the original network. This approximation becomes exact in the limit of 
networks with an inhnite number of neurons, the so called thermodynamic limit. Clearly this means that 
for hnite-size networks, the mean-held theory provides only an approximation of the real behavior of the 
system, and therefore may neglect important phenomena. 

for a class of random neural networks it was proved that in the thermodynamic 


For example, in 16 


limit the system undergoes a sharp transition from a static regime to chaos, while in the hnite-size case 
this occurs gradually by a cascade of transitions that are not displayed by the mean-held approximation. 
Another important example is shown in 17 , where the authors proved that a stochastic hnite-size network 
can exhibit arbitrarily high levels of cross-correlations of the neural activity, while this phenomenon 
disappears when the size of the system grows to inhnity. Clearly, these macroscopic differences in the 
dynamical and statistical behavior of hnite and inhnite-size networks may have important consequences 
on the information processing capability of the system. In other terms, the oversimplihcation of the mean- 
held approximation may hide important neural processes that are fundamental for the comprehension of 
the brain. 

This explains the number of mathematical techniques that have been developed to quantify hnite-size 
ehects beyond the mean-held theory, such as the linear noise approximation 18 , the density functional 
approach 19 , large-deviations theory 20 , path-integral methods 21 , etc. Typically, these hnite-size 


methods can be applied to networks composed of a hnite but large number of neurons. However, neural 
masses in a cortical minicolumn contain only few tens of neurons since, according to different estimates 
22 - 24 , a minicolumn in primates contains about 80 — 100 neurons. Thus, for such small networks the 


mean-held description and the previously introduced techniques turn out to be inappropriate, since their 
hnite-size effects can be relevant. 

The analysis of the dynamics of small neural networks was started by Beer 


25 , who studied the 


bifurcations of networks of arbitrary size in highly symmetric cases, namely with rigid constraints on 
the strength of the synaptic weights. Through this analysis, he showed that small networks can exhibit 
complex dynamical behavior, which may have important neurobiological implications. In our article we 
extend his analysis to a more biologically plausible network of arbitrary size without specihc conditions 
on the synaptic weights, and in particular we show another important difference between a finite-size 
network and its mean-held counterpart. In more detail, here we consider a deterministic hnite-size hring- 
rate network model with excitatory and inhibitory populations composed of arbitrarily few neurons (we 
consider Ne neurons in the excitatory population and Ni neurons in the inhibitory one), which are 
indistinguishable within each population. Then, we perform a numerical analysis of the dynamics that 
emerge by varying the external input current to the network and the strength of inhibition, and we 
introduce a mathematical theory that allows us to describe local bifurcations analytically. In particular, 
we obtained macroscopic differences with the mean-held approximation when the system has strong 
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inhibitory synaptic weights. In this case, through a phenomenon of spontaneous symmetry-breaking, 
the neural network undergoes a special bifurcation known as branching point 26 , from which multiple 


solutions of the neural equations emerge. On the new branches, new bifurcations can occur, enriching 
considerably the complexity of the bifurcation diagram of the neural network. For this reason, symmetry 
and branching points bifurcations may reveal important mechanisms through which neurons interact and 
process information, which are not accounted for by the mean-field approximation. Therefore a detailed 
analysis of their properties may be very important in order to understand the role of the mesoscopic scale 
in the emergence of brain’s highest functions. 

It is important to observe that since we consider identical neurons within each population, our model 
lies in the context of the bifurcation theory of dynamical systems with symmetry, which is known as 
equivariant bifurcation theory 27 . This usually requires the reader to be familiar with advanced notions 
of group theory. Notwithstanding, here we follow a simpler approach, so the mathematical analysis 
introduced in this article is self-contained and accessible also to non-technical readers. 

Symmetry-breaking and branching point bifurcations have already been studied in neuroscience 
through equivariant bifurcation theory, but in a conceptually different way. For example, in the the¬ 
ory of visual hallucinations developed by Ermentrout and Cowan 28 , symmetry was used to evaluate 


hallucinations in primary visual cortex under the effect of drugs. They idealized the cortex as a plane 
and described the local activity of a population of neurons by means of neural field equations, so that 
each population is univocally identified by its position in the continuous space of the plane. Their theory 
exploits the symmetry the neural field equations inherit from the geometrically regular structure of the 
anatomical connections in primary visual cortex. In equivariant bifurcation theory, this symmetry is 
described by the invariance of the equations under the action of the Euclidean group E (2), which is the 
group of rigid motions in the plane, generated by translations, rotations, and reflections. 

However, the network analyzed in our work is made of two populations containing a finite number of 
neurons, which are identified by a discrete index. In each population the neurons are identical, therefore 
the equations are symmetric under permutations of the neural indexes within a population. In more 
mathematical terms, our equations are invariant under the action of the group Snj^. x where Sn^ is 
the permutation group on Na items (also known as symmetric group). This is a completely different kind 
of symmetry compared to that of the Euclidean group E (2), and it allows us to study in an analytically 
simple way networks made of a finite number of neurons. Indeed, this is conceptually different from 
the theory of Ermentrout and Cowan, which relies on infinite-dimensional neural field equations. So 
their theory does not describe finite-size effects, and does not use the underlying symmetry for this 
purpose. In this respect, our theory is more similar to that developed in 29-32 , where the authors 


exploited the invariance under the symmetric group in a system with all-to-all coupling, and the 
subsequent spontaneous symmetry-breaking, as a possible mechanism to explain the origin of animal 
species. However, we are not aware of any application of this symmetry to the study of bifurcations in 
neural networks described by specific equations, such as Hopfield 33 or Wilson-Cowan 34 ones. 

This article is organized as follows. In Sec. ([^ we describe the neural model we use. Then we 
start Sec. by explaining intuitively the main topic of this article, namely the formation of new 
branches of solutions through spontaneous symmetry-breaking, depending on the strength of inhibition 
(see SubSec. ( |3.1[ )). This is followed by a numerical and analytical study of the bifurcation points of the 
network in weak and strong-inhibition regimes (see SubSecs. (3.2) and (3.3) respectively). As we will 
see, the complexity of dynamics depends also on the size of the inhibitory population, therefore we start 
our analysis by showing the behavior of the network in the simplest case, namely Nj = 2, and then we 
explain how to generalize these results to a larger inhibitory population (see SubSec. (3.4)). We also 
show why the mean-field theory is an oversimplified description of the network (SubSec. (3.5)), and that 
the finite size effects studied in this work for a fully-connected network may be even stronger for more 
realistic topologies of the synaptic connections (SubSec. (3.6)). In Sec. (|^ we examine the importance 
and the biological implications of our results. Finally, more details about the analytical calculations are 
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shown in Supplementary Materials. 


2 Materials and Methods 


The analysis of the dynamics of a neural network is a mathematically hard problem. Its complexity 
depends on the degree of biological plausibility of the system, which in turn is determined by: 

• the variety of neural populations and realistic ratios between their sizes; 

• the finite size of the network; 

• the random behavior of the system; 

• a realistic topology of the synaptic connections; 

• plausible single-cell models for the soma and the synapse; 

• delays in the transmission of the electric signals through the axons. 

Among the points listed above, we select several biological plausible features that lead our model to be 
characterized by: 

• two neural populations of excitatory and inhibitory neurons; 

• a generic finite number of neurons in the network and in each population; 

• non-symmetric and arbitrarily strong synaptic weights. 

For the sake of clarity, in this article we perform a numerical and analytical analysis in the case of 
just two neural populations, since they are sufficient for describing dynamical behaviors, such as neural 
oscillations, which are thought to be at the base of fundamental cognitive processes. It is important to 
observe that our analysis can be extended to a generic number of neural populations, as we discuss briefly 
in SubSec. (S4.1) of the Supplementary Materials, but this problem will be addressed in more detail in 
another article. In order to make the analysis of the two-populations network analytically tractable, we 
make some simplifying assumptions that considerably reduce the complexity of the problem, without 
losing the most fundamental properties of a neural network: 

• non-spiking neurons; 

• identical neurons in each population; 

• dense connections (i.e. fully-connected topology); 

• absence of axonic delays; 

• absence of random noise. 


These features are supposed to not compromise the plausibility of the final result. Actually, this work is 
intended as a model of a neural mass within a cortical column, where neurons are homogeneous and 
the density of the connections is known to be high, according to anatomical data (see for example [^). 
This justifies the use of identical neurons in each populations and of a fully-connected topology. Moreover, 
since a single column is localized in a mesoscopic area of the brain, delays can be neglected. However, 
following 35,36 , our analysis can be extended to the case of delay differential equations, if desired. Also 
random noise may be taken into account, as shown in 17 , but this will be studied in detail in another 


article. To conclude, the choice of spiking or rate models is still under debate in the scientific community 
(see for example 37-^) so this, together with the need for a mathematically tractable model, justifies 
our decision to work with a rate neural network. 

In more detail, we consider a widely used model of rate network [^[T^[T^[T7|[^[3^ : 
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where N > 2 represents the number of neurons in the network. The function Vi (t) is the membrane 
potential of the ith neuron, while Ti is its membrane time constant. Mi, the total number of 
connections to the ith neuron, is a normalization factor that has been introduced to avoid the explosion 
of the total synaptic input i^j (0) equation Q in the thermodynamic limit N —>■ oo. 

The matrix J, also known as structural or anatomical connectivity, represents the specification of all the 
synaptic wirings that are physically present between neurons. It also quantifies the strength of these 
connections. So is the synaptic weight from the jth to the ith neuron, and for simplicity it is 
supposed to be deterministic and constant in time, for every pair of neurons. £/j (•) is the activation 
function of the jth neuron and converts its membrane potential into the corresponding firing rate 
i/j = £/ (Vj). S'-shaped (i.e. sigmoidal) activation functions are biologically plausible, however usually 
piecewise-linear functions are used in order to find analytical results 41-43 . In our work we show how 


it is possible to obtain analytical expressions for the equilibrium points and the local bifurcations of the 
network when M (V) is the so called algebraic activation function: 


M (V) 


1 + 


¥iv-vn 


( 2 ) 


Here is the maximum firing rate, A determines the slope of the activation function when 
is fixed, and is the horizontal shift, so it can be interpreted as a firing threshold for the membrane 
potentials. In Eq. 0 . li are external input currents, which in this article are supposed to be constant in 
time, in order to perform the bifurcation analysis of the network. 

In order to make our analysis analytically tractable, we suppose that all the parameters of the system 
are indexed only at the population level. In other terms, within a given population, or between two 
fixed populations in the case of the synaptic weights, the parameters are homogeneous. If we define Ne 
(N j) to be the size of the excitatory (inhibitory) population, and if we suppose that the neurons with 
indexes i = 0,..., Ne — 1 belong to the excitatory population, and those with i = Ne, ..., N — 1 (with 
N = Ne + Nj) to the inhibitory one, this means that the synaptic connectivity matrix is structured as 
follows: 


^EE ZeI 1 
ZiE Zll J ’ 




Jaa (IjVc — IdjV„ ) , 

Jay I Vq ,N a , 


for a = P 
for a ^ d 


where the block Zap is a Na x Ng matrix that represents the connections from the population P to the 
population a, with a, P € {E,I}. Moreover, is the N^ x Ny all-ones matrix (here we use the 


def • 


simplified notation '= while is the Na x Na identity matrix. Since, according to 

experimental mea surements, about 80% of neocortical neurons are excitatory and the remaining 20% 
are inhibitory 44 , in this article we consider ^ = 4. The zeros on the diagonal lines of the matrices 


Zee and Zii are due to the absence of self-connections in biological networks. The real numbers Jee, 
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Figure 1: Structure of the network considered in this article. 


Jiij Jei, Jie are free parameters that describe the strength of the interactions between and within the 
neural populations. Clearly we have Jee, Jie > 0, and J//, Jei < 0, which also means that 
Me = Ml = N — 1. It is important to observe that, compared to 


15 , we have chosen a different 


normalization of the synaptic weights. In both the cases, the total synaptic input Jij'^j (Yj (^)) i® 

convergent in the thermodynamic limit, but they behave in different ways when the network is mixed 
with both finite-size and infinite-size populations. It is easy to check that according to our 
normalization, if one population has inhnite size, then the contribution of the finite-size population in 
the total synaptic input goes to zero, which seems reasonable to us. On the other side, with the 
normalization introduced in 


15 , the infinite-size and finite-size populations both provide finite 


contributions to the total synaptic input. However, if required, it is possible to switch from our 
normalization to the other one by performing the following substitution: 


JaP —I -TT- JaP 

in all the analytical results that we will obtain in the next sections and in the Supplementary Materials. 

Moreover, from our assumption on the indexes, we obtain that the external input currents are divided 
in two vectors, Ie and Ij, namely the inputs to the excitatory and inhibitory populations. Clearly we 
get: 


Ic 


lalNc 


where is the x 1 all-ones vector. Therefore the structure of the network is that shown in Fig. 0. 
The same subdivision between excitatory and inhibitory populations is performed for the other parameters 
of the network, namely r, A, Vr- 

To conclude, for the sake of clarity we observe that under our assumption on the neural indexes, we 
can rewrite the system Q in the following explicit way: 


E : 


= + M i = o....,WE-i 


3=0 




iVR, - 1 




(3) 
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Population Size 

Synaptic Weights 

Activation Functions 

Other 

00 

II 

o 

II 

max max -i 

— -L 

Te =TI = 1 

Ni = 2 

Jei = —70 

Ae = A/ = 2 



Jie = 70 

Vj =Vf ^2 



Table 1: Values of the parameters used in this article. 


3 Results 

The bifurcation analysis we perforin provides an overview of the dynamics the model is able to exhibit, 
depending on two parameters: the static external currents Iej in (|^. In particular, we present this 
analysis for increasing values of a third parameter: the self-inhibition strength J//. We focus on this 
parameter, instead of the other synaptic weights (i.e. Jee, Jei, Jie), since J// is directly responsible for 
the main phenomenon analyzed in this article, namely the formation of the branching point bifurcations. 
For the remaining synaptic weights we will provide only a qualitative description of the effects they exert 
on the system. The non-varying network parameters are chosen as in Tab. 0. In particular, we consider 
a network made of Ne = 8 excitatory neurons and Nj = 2 inhibitory ones, since we want to study 
finite-size effects in small neural masses. 

We perform a detailed bifurcation analysis by means of numerical tools and, when possible, through 
analytical techniques. The numerical analysis is performed by exploiting the CLMatCont Matlab toolbox 
45 and XPPAUT 46 , that are grounded in the mathematical theory of bifurcations described in 26p7 , 
while the analytical results are based on elementary methods from linear algebra. It is important to 
underline that bifurcations are defined by many conditions. Nonetheless, in our analytical study we 
checked only the conditions on the eigenvalues of the network, since they proved sufficient to reproduce 
the numerical results. Due to the high variety of the bifurcations the system exhibits, a rigorous check 
of all the remaining conditions is beyond the purpose of this article, and is left to the most technical 
readers. 


3.1 Intuitive interpretation of the branching points 


In mathematics 
theory 


27 


the branching point bifurcations are described by the so-called equivariant bifurcation 
namely the study of bifurcations in symmetric systems. Being the latter rather technical, 
here we prefer to follow a more intuitive approach to the problem. So, first of all, we have to observe 
that according to bifurcation theory, local bifurcations are calculated by means of the eigenvalues of 
the Jacobian matrix of the network, evaluated at the equilibrium points. Therefore, we start by setting 
= 0 Vi in Eq. The main observation of this article is that the system shows two different 
behaviors depending on the strength of J//. As far as inhibition is weak (this will be quantified rigorously 
below), the equilibrium points in each population are homogeneous: 


fJ- = 


f times Nj—times \ 

fJiE ) • ■ • ; f-^E 7 /^/ 5 • ■ ■ 7 


(4) 
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where fiE and are the solutions of the following system of algebraic non-linear equations, obtained 
from 


( ^ Jee-b^e (pe) + -^£iJeis^i {pi) + Ie = Q 

\ (5) 

{pE, pi) + 'W^’^IE-^E {pe) + JiiS^I (pi) + Ii = 0 

The curves defined by Eqs. ^ {x, y) = 0 and ^ {x,y) = 0 V (x, y) e are the so called nulldines of 
the network. Fig. ([^ (top) shows an example obtained for J/j = —10, while the remaining parameters 
are chosen as in Tab. Q. In Sec. (S2) of the Supplementary Materials we show how to get approximate 
analytical solutions for pej- 

From (§, we can see that the Jacobian matrix JT of the network on the primary branch of solutions 
Q is: 


^=[ 


r o^ee 

Jei ] 

a ^ _ J 

[ JlE 

,/ri J’ 



- -I- {pa) (liVa - Idiv„), for a = /3 


( 6 ) 


{pft) lN^,Nfi , for a 7^ d 

therefore we can prove (see SubSec. (S3.1) of the Supplementary Materials) that its eigenvalues are: 


5 + ri± \j{5 — v)^ + 47 
Ao,i = -r-, Xe 




te 


Xi = — 


1 Jii 
Vi ^ N-1 


{pi) 


(7) 


where: 


5 = -h — :^Jees/e {pe) , v = -h — :rJiii^i {pi) , 7 = . V '^eiJieJ^^e {pe) {pi) (8) 

According to bifurcation theory, the system undergoes special bifurcations when one of its eigenvalues is 
equal to zero. In particular, the branching point bifurcations are given by the condition A/ = 0, so this 
allows us to define the weak and strong-inhibition regimes quantitatively by the conditions A/ < 0 and 
Xi > 0, respectively. 

When the network undergoes a branching point bifurcation, we observe the formation of heterogeneous 
membrane potentials in the inhibitory population: in other terms, under strong inhibition the symmetry 
of the system is broken. Intuitively, this can be understood as in Fig. When |J//| is small, there 
is only one valley or basin in the “energy landscape” of the network. On the other side, for strong 
inhibition we observe the formation of multiple valleys, and a small perturbation determines to which 
one the inhibitory potential will converge. For this reason, now multiple new branches of the equilibrium 
points emerge, which are described by the following stationary solutions: 


—times 


p = 


Pe, •.., Pe, pi,o, ■ ■ ■, Pi,Nj-i 
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Figure 2: The two figures on the top show the nullclines of the network for fixed values of Ie,i in a weak- 
inhibition regime, obtained for Jn = —10 and the values of the parameters in Tab. Q. Their intersection points 
(black dots) correspond to the solutions of the system §. The top-left figure was obtained for Ie = 10, // = —10, 
while the top-right figure for Ie = 13, 7/ = —10. Clearly the system § admits multiple solutions for specific 
values of the parameters. The two figures at the bottom show the solutions jib and /r/ (bottom-left and bottom- 
right panel, respectively) of the system nt for the same values of the parameters, but with varying current Ie- 
The black curve represents the primary branch of the network equations, and for Ie = 10 and 7 _b = 13 it admits 
one and three solutions respectively (see the black dots at the intersection with the vertical dashed lines). The 
{jj,E,iJ,i) coordinates of these solutions correspond to those of the black dots in the top panels of the figure. Here 
we do not care about the stability of the solutions, which is shown for the sake of completeness in Fig. (S5) of 
the Supplementary Materials. 
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A/ < 0 


A/ = 0 


A/ > 0 





Figure 3: Spontaneous symmetry-breaking. For weak inhibition, the eigenvalue Aj (see Eq. 0) is negative, 
and the system (|^ has only one stationary solution. This solution is stable and symmetric, but for increasing 
inhibition A/ changes sign, therefore the solution becomes unstable and the network chooses randomly between 
two new alternative stable states, breaking the symmetry. This phenomenon can be understood intuitively by 
means of a ball that rolls in a double well potential in order to reach a state of minimum energy, that we interpret 
as a stable stationary solution of Eq. 0- 


where and are the solutions of the following system of algebraic non-linear equations, obtained 
from 0 in the stationary regime: 


Nj-l 

JeeJ^e (ij-e) + ^ M (pij) -f /e = 0 


— + ■^3[Jies^e (me) -|- = 0 for i = 0,..., A^/ — 1 

j=o 


( 10 ) 


For example, for Nj = 2 Eq. (101 can be written more explicitly as follows: 


^ (me, Mr.o, M/,i) — y^ME + Jee-^e (me) + (m/.o) + (m/.i)] + /e = 0 

< {^E, + ^3[JieS^E {iJ-e) + ~ ^ ( 11 ) 

M’ (me, Mr. 0, Mr,i) "== —y^Mr.i + ^-eiJieJ^e (mb) + (Mr.o) -f // = 0 


The surfaces {x,y,z) = 0, {x,y,z) = 0, J^{x,y,z) = 0 'i{x,y,z) G are a higher-dimensional 

extension of the nullclines ^ {x, y) = 0 and (a;, y) = 0 V {x, y) G that we encountered in the weak- 
inhibition regime. Sometimes they are called nullsurfaces (see for example [25| ). Fig. Q (top) shows an 
example obtained for J// = —100, while the remaining parameters are chosen as in Tab. Q. 

For the sake of clarity, here we treat in detail only the case Nj = 2, while we will show some results 
for Nj > 


2 in SubSec. (3.4). 
secondary branches 0 is: 


So, for Nj = 2, from Eq. (11) we get that the Jacobian matrix on the 
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Figure 4: The figure on the top shows the nullsurfaces of the network for fixed values of Ie,i in a strong- 
inhibition regime, obtained for Jn = —100, Ie = 5, 1/ = —10 and the values of the parameters in Tab. 0- The 
black intersection point corresponds to the solution of the system ( |11[ | on the primary branch, while the purple 
dots represent the solutions on the secondary branches. The two figures at the bottom show the solutions fiE 
and /ij (left and right panel, respectively) of the system (111 for the same values of the parameters, but with 
varying current Ie- The black and violet curves represent respectively the primary and secondary branches of the 
network equations. For Ie = 5 the system admits three solutions {ijle,IIi) (see the dots at the intersection with 
the vertical dashed fines: we have just three solutions because /rj has three intersection points, two of which, i.e. 
the pink ones in the right panel, correspond to the same /i_B, i-e. the pink dot in the left panel). The {fiE,fJ‘i) 
coordinates of these solutions correspond to those of the dots in the top panel of the figure. Again, here we do not 
care about the stability of the solutions, which is shown for the sake of completeness in Fig. 0 and in Fig. (S9) 
of the Supplementary Materials. 
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^10 = (mb) Ibe , ^11 = —, ^12 = ^ -e^/ (m/.i) : 


^20 = {( J - e ) 1 a ) e , ^21 = ( M /. o ) , ^22 = — 


Tl 


(here T is the transpose of a matrix), as explained in more detail in SubSec. (S4.1) of the Supplementary 
Materials. Intuitively, on the primary branch the Jacobian matrix was a 2 x 2 block matrix (see Eq. (®), 
because we had only an excitatory and an inhibitory membrane potential {fj-E and /r/ respectively), while 
on the secondary branch it is a 3 x 3 block matrix, because now the two inhibitory neurons have different 
potentials (/i/.o and The Jacobian matrix on the new branches can be calculated for a network 

with a generic number of inhibitory neurons (see the Supplementary Materials for more details). 

Finally, we observe that it is possible to find relations between the inhibitory membrane potentials in 
the strong-inhibition regime, which prove very useful when we calculate analytically the local bifurcations 
of the system. For example, in the case Nj = 2, from the second and third equation of the system 0- 
after some algebra we get the following forth-order polynomial equation: 




(13) 


whose coefficients depend on as follows: 
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Eq. (131 can be solved analytically, providing an explicit expression of as a function of ^/, 0 ) which 
will be used in SubSec. (3.3) to evaluate the local bifurcations in the strong-inhibition regime. 


3.2 Weak-inhibition regime (A/ < 0) 

As we said, we want to understand how the network’s dynamics changes when we vary the external input 
currents Ie.i and the strength of the synaptic weight J//. For the sake of clarity, in this section we show 
the results that we obtain when we vary these parameters one by one, because this allows us to introduce 
the concepts of codimension one and codimension two bifurcation diagrams. 

We start by considering fixed Ie,i and inhibitory strength J// = —10. As we can see from Fig. ([^ 
(top), Eq. © admits multiple solutions, depending on the specific value that we have chosen for the 
current Ie- Then, we start to vary Ie continuously, while keeping Ij and J// fixed. In this way we 
can plot how the solutions of the system ([^ change as a function of I^;. As we anticipated in 
SubSec. (3.1), the curve that we obtain is called the primary branch of the network, see Fig. ([^ (bottom). 
This hgure represents the codimension one bifurcation diagram of the network, from which we see there 
are two special points, called local saddle-node bifurcations, or LP for short. They identify the value of 
Ie for which the number of solutions of Eq. (§ changes (notice the correspondence with the top panels 
of Fig. (§). These are the first example of (local) bifurcation points that the neural network exhibits, 
and that lead to the formation of hysteresis (see the left panel in Fig. (©)• Hysteresis was suggested to 
describe short-term memory, since a sufficiently strong input can lead the system to reach a stable-high 
level of activity that is maintained when the input is turned off 48 -50 . This phenomenon, known as 


reverberation, namely the persistence of neural activity sustained internally in the brain after a stimulus 
is removed, can be achieved through bistability, which can be present intrinsically at the single cell level 
or generated by recurrent excitatory connections. A typical example is the working memory, which is 
intimately related to the prefrontal cortex and that can retain information for a time span of the order 
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Figure 5: 
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integrator 
relevance. 


On the left, we show an example of hysteresis displayed by the system. The plain lines describe stable 
while the dashed line the unstable ones. On the right, we show an example of catastrophe manifold, 
highlights three different behaviors of the network for increasing values of |//|: leaky integrator, perfect 
and switch (red curves). The readers are referred to [M 55 for more details about their biological 


of seconds, before it is erased by events such as noise or distracting stimuli. It is important to observe 
that even if reverberation requires bistability, the latter can be present without hysteresis, but very often 
they coexist, as in our model. 

If now we continuously vary both Ie and //, while keeping J// fixed, we can plot the solutions of Eq. (§ 
as a function of both the external currents, so now ^e,i = ^^E,I [Ie-: Ii) define three dimensional manifolds 
(see Fig. ([^ for ^le)- On these manifolds, the special points LP depend also on //, therefore they form 
a set of points called saddle-node curves (see the blue curves in Fig. (©)• For visual convenience, these 
curves are projected on the Ie — Ii plane (see Fig. 0). defining the so-called codimension two bifurcation 
diagram of the network (which is, from the analytical point of view, our main interest in this article, 
while the codimension one diagram is partially calculated in Sec. (S2) of the Supplementary Materials). 
However, from Fig. 0 and the bottom panels of Fig. Q, we can see that the saddle-nodes are not the 
only bifurcations we have in our network. For some values of the pair Ie — Ii, the system can undergo 
also a so called local Andronov-Hopf bifurcation, or H for short. These bifurcations are represented by 
the red curves in Figs. ([^ and ([^, and correspond to the emergence of neural oscillations, which are 


thought to play a key role in many cognitive processes 51 


Hereafter, we list all the bifurcations the system undergoes, dividing them in groups depending on 
the codimension of the bifurcation, which is defined as the number of parameters (in our case Ie,i) that 
must be varied for the bifurcation to occur. Although only few of them are represented in Fig. ([^, the 
complete set of codimension one bifurcations our system undergoes is: 

• local saddle-node bifurcations (LP), for which new equilibria arise, or collide and annihilate each 
other; 


local Andronov-Hopf bifurcations (H), where stable or unstable self-sustained oscillations, described 
by limit cycles, arise or disappear; 

global homoclinic bifurcations, where limit cycles vanish in a particular equilibrium point (i.e. a 


neutral saddle, see SubSec. (3.3|, or a saddle-node), giving rise to an orbit with infinite period; 


global limit point of cycles bifurcations, at which new limit cycles arise, or collide and annihilate 
each other. 
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Figure 6: Three-dimensional bifurcation diagram equivalent to that in Fig. 0. The values of are shown 
on the z—axis as functions of Ie — Ii- Here we plot only the local bifurcations (blue: saddle-node curves, red: 
Andronov-Hopf curves) that bound the gray/red regions representing the stable/unstable equilibrium point areas. 
All the codimension two bifurcations are reported in Fig. Q. 


The codimension one diagrams collecting all these bifurcations are shown in Fig. (S5) of the Supplemen¬ 
tary Materials. Moreover, on the curves defined by these bifurcations, and that are obtained by varying 
both Ie and //, the following codimension two bifurcations appear (see Fig. [^: 

• local cusp bifurcations (CP), on the saddle-node curves; 

• local generalized Hopf bifurcations (GH), that divide subcritical Andronov-Hopf curves from super¬ 
critical ones; 

• local Bogdanov-Takens bifurcations (BT), that represent the contact point between the saddle-node, 
Andronov-Hopf (ending here) and homoclinic curves; 

• global cusp of limit point of cycles (CPC), on the limit point of cycles curves; 

• global saddle-node on invariant circle (SNIC), where a saddle-node bifurcation occurs simultaneously 
with a homoclinic bifurcation. 

It is worth remarking that the bifurcation diagram in Fig. ([^, obtained from the voltage-based model 
0 in the weak-inhibition regime, is qualitatively similar to that of the activity-based Wilson-Cowan 


model (see Fig. 2.12 in 52 ). These two kinds of models are obtained from neural mass equations 


through two slightly different hypothesis about the postsynaptic potentials 0. For this reason, in the 
literature it has always been assumed implicitly that the two models can exhibit qualitatively similar 
dynamics. The strong resemblance of their codimension two bifurcation diagrams proved in this article 
confirms rigorously this intuition (notwithstanding, in the next section we will show that things may 
change significantly in the strong-inhibition regime if we take into account finite-size effects). 

Interestingly, Fig. presents two of the so-called catastrophe manifolds [53] , one of which is shown 
in the right panel of Fig. (0. This figure emphasizes the ability of the model to describe three different 
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Figure 7: Complete codimension two bifurcation diagram on the Ie — Ii plane in the weak-inhibition regime. The 
Andronov-Hopf bifurcation curves (red lines) are divided into supercritical (plain) and subcritical (dashed) por¬ 
tions. The supercritical/subcritical portions are bounded by a Generalized Hopf bifurcation, GH, and Bogdanov- 
Takens bifurcations, BT. The latter are the contact points among saddle-node bifurcation curves (blue lines), 
Andronov-Hopf bifurcation curves (red lines), and homoclinic bifurcation curves (hyperbolic-saddle/saddle-node 
homoclinic bifurcations are described by plain/dashed orange curves). SNIC bifurcations identify the contact 
point between the saddle-node curve and the homoclinic one. From GH originate two Limit Point of Cycles 
curves (dark green lines) that collapse into the homoclinic curves. Before this, they present a cusp bifurcation, 
CLC. Each saddle-node curve shows, in addiction to BT, a cusp bifurcation, CP. 
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behaviors: leaky integrator, perfect integrator and switch. This triad represents the main ingredient 


for describing a mechanism which was proposed to explain interval timing by neural integration 54 55 


According to the theory, by changing some parameters of the network, it is possible to modify the shape 
of the equilibrium curve on the catastrophe manifold, generating a ramping activity that can explain 
Weber’s law of time perception 56 . This phenomenon can easily occur in our model, where the shape 


of the equilibrium curve can be changed dynamically by varying the input currents. 

Now we want to prove some of our previous results on bifurcations from the analytical point of 
view. Often, in dynamical systems, global bifurcations are harder to study analytically, therefore here 
we focus on LP, H and BT. So first of all we observe that, according to bifurcation theory 26,47 , LP is 


mathematically defined by the condition that one real eigenvalue of the Jacobian matrix becomes equal to 
zero. Therefore, from Eq. Q, we conclude that for our network this bifurcation occurs whenever Aq = 0 
or Ai = 0, because Xe is always negative, while A/ = 0 defines the branching point bifurcation. For 
example, if <5 + ry < 0, the condition Aq = 0 is equivalent to 6r] = j which, according to Eq. ® , provides: 
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where we have defined the parameter D = ^e- Now we invert ip-i) (more details are provided in 
SubSec. (S3.2.1) of the Supplementary Materials), obtaining: 


and from Eq. ([^ we get: 
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and they define analytically the blue curves in Fig. Q. As we said, this is not sufficient to prove 
that Eqs. (14) - ( [T^ describe saddle-node curves, since we should check also the corresponding non¬ 
degeneracy conditions. Nevertheless, we observe a perfect agreement between these analytical curves and 
those obtained numerically by CLMatCont and XPPAUT, therefore for simplicity we do not check the 
remaining conditions and we leave them to the most technical readers. We adopt the same approach for 
the remaining bifurcations we are about to describe. 

they appear whenever the network has 


Now we focus on the H bifurcations. According to 26,47 


conjugate purely imaginary eigenvalues. Since Xe,i are always real, this condition can be satisfied only 
by Aop, by setting J -|- 77 = 0 and (<5 — -I- dy < 0. In particular, from the equation S + ij = 0 we get: 
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where D = as before. Following the same procedure introduced before for the LP curves, we obtain 
a set of parametric equations for the pairs Ie — 11 that generate the H curves, with parameter o € 
[0/,Dd] U [Dc,0e], where: 
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Now, as we said, BT represents the point where the LP and H curves meet each other, and identifies 
also the end of the H curve. Clearly, from the condition Aq = 0 or Ai = 0 that defines the LP curves, 
and the condition Ao,i = itw (where i represents the imaginary unit) that defines the H curves, we get 
Aq = Ai = 0. This is the condition that defines analytically the BT points, or equivalently Dbt = Oe,/ 
as given by Eq. ( [T^ , from which the coordinates of the BT points in the Ie — Ii plane can be easily 
obtained through Eq. (§. The remaining local bifurcations (i.e. CP and GH) are analytically intractable, 
therefore we cannot study them beyond the numerical results. 

To conclude this subsection, now we describe briefly the effect of the variation of the remaining 
synaptic weights on the codimension two bifurcation diagram, considering the weights in Tab. Q as 
reference point. As we said, their variation does not generate interesting phenomena such as the branching 
point bifurcation which is obtained by varying J//, so we dedicate less space to these parameters. 

For Jee ^ 10, the two LP curves become larger and larger on the Ie axis (i.e. the distance between 
their vertical asymptotes increases). Moreover, the curves get closer and closer to each other, by shifting 
on the Ie axis, until they intersect and their oblique parts (i.e. those between the BT points) overlap. If 
we increase Jee further, the LP curves split again in two disjoint parts, each one presenting two BT and 
two CP bifurcations (so the total number of CP points increases from two to four). Between each pair of 
BT points (on the same LP curve) there is a H curve. These curves are very close to the corresponding 
LP curves, and if we increase Jee further they disappear, together with the BT bifurcations. So for 
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very large Jee we get only two disconnected LP curves, or in other terms for very strong excitation the 
oscillatory activity cannot be sustained anymore. Also on the opposite side, namely for weak inhibition 
(i.e. Jee —>■ 0), the H curves disappear. Moreover, the width on the Ie axis of the LP curves decreases, 
i.e. the distance between their vertical asymptotes becomes smaller and smaller, until the asymptotes 
collapse on each other and the LP curves disappear as well. 

For \Jei\ ^ 70, the width of the two LP curves remains almost constant, while the distance between 
them (and therefore also the length of the H curves) increases continuously. On the other side, for 
\Jei\ 0, the two LP curves get closer and closer to each other, until they intersect and then their 
oblique parts between the BT points overlap. If we decrease \Jei \ further, similarly to the case with large 
Jee, the two LP curves split in two disjoint parts, each one presenting two BT and two CP bifurcations. 
For even smaller values of \Jei\, the BT, CP and H bifurcations disappear, while the LP curves disappear 
for \Jei\ = 0. 

For JjE ^ 70, the LP curves are stretched vertically and shifted downwards along the // axis. Clearly, 
in the opposite direction (i.e. for Jje —> 0), they are compressed, and while the Ie coordinates of the 
vertical asymptotes remain almost unchanged with Jie, the two CP points get closer and closer to each 
other. At the same time, the two H curves tend to overlap. At some value of Jje, the two CP bifurcations 
touch each other and disappear, so the two LP curves become tangent. If we further decrease Jje, the 
two LP curves split again, one over the other, and the BT and H bifurcations disappear. 

All the phenomena that we have just described are qualitatively similar for different values of J//, so 
they occur also in the strong-inhibition regime for the primary branch. 


3.3 Strong-inhibition regime (A/ > 0) 


In the strong-inhibition regime (in particular here we consider the cases Jjj = —34 and J// = —100), 
most of the features of the weak-inhibition bifurcation diagram are preserved. However, besides the 
bifurcations explained in SubSec. ( |3.2[ ), from Figs. Q, ([^, (10), we can see that the system undergoes 
also the following codimension one bifurcations: 


• local branching point bifurcations (BP), at which two or more equilibrium curves intersect each 
other; 


• local torus bifurcations (TR), at which the limit cycles are characterized by a quasi-periodic motion; 
and the following codimension two bifurcations: 

• local zero-Hopf (neutral saddle) bifurcations (ZH), that involves the Andronov-Hopf curves and, 
in our case, the branching point curves. Around this point, both subcritical and supercritical 
Andronov-Hopf bifurcations exist. 


In particular, the branching point bifurcations lead our model to show multiple branches of stationary 
solutions for suitable current values. This is a finite-size effect, due to the finite number of neurons in 
each population, that leads to a richer set of Jacobian matrix eigenvalues than that obtained by using 
methods based on the reduction of the number of equations, such as the mean-held approximation (see 
SubSec. (3.5) for more details). In order to thoroughly investigate the bifurcations the system undergoes 
in presence of strong inhibition, we start by analyzing the codimension one bifurcation diagram for 
Jii — —34. In particular, the diagram in Fig. § is obtained by letting // = —10. It turns out that, 
in addition to the primary equilibrium point curve (black line), new branches of stationary states (violet 
lines) emanate and collapse in two branching point bifurcations, BP. These secondary branches hold 
supercritical Andronov-Hopf bifurcations that give rise to stable limit cycles. Instead, on the primary 
branch, we hnd two saddle-node bifurcations and a subcritical Andronov-Hopf bifurcation, whose unstable 
limit cycles vanish in a homoclinic orbit. We remind that branching point bifurcations occur because 
A/ changes sign. We also observe that for Nj = 2 the inhibitory neurons have the same bifurcation 
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Figure 8: Codimension one bifurcation diagrams for (left) and /r/ (right) as a function of Ie, for Ju = —34 
and Ii = —10. In both the graphs, the stable/unstable primary equilibrium curve is described by plain/dashed 
black curves, while the secondary ones are described by plain/dashed violet curves. Andronov-Hopf bifurcation 
appear on both the primary and secondary equilibrium curves, giving rise to unstable and stable limit cycles 
respectively. In particular, the unstable cycles collapse into a homoclinic bifurcation, described by an orange line. 


diagram (see Fig. ([^, right). However, this does not mean that the inhibitory membrane potentials are 
homogeneous. Indeed, when an inhibitory neuron is on the upper secondary branch (see the violet curve 
above the primary equilibrium point curve in Fig. (|^, right), the other one is in the lower secondary 
branch, so indeed they are heterogeneous. 

For Jjj = —100, secondary branches of equilibrium points are still present, see Fig. (§. Together 
with a subcritical Andronov-Hopf bifurcation, they unveil also saddle-node ones. In particular, the former 
generates unstable limit cycles that become stable after having crossed the limit point of cycles bifurcation 
(dark green line). For increasing values oi Ie, the stable limit cycles collapse into the unstable limit cycles 
originated from the subcritical Andronov-Hopf bifurcation belonging to the primary equilibrium point 
curve. Before collapsing, the stable limit cycles undergo torus bifurcations (gray line). 

By varying also //, we obtain the codimension two bifurcation diagrams displayed in Fig. ( [l0| . It 
is worth noting that the branching points the system undergoes generate two bifurcation curves (light 
green lines) that pass through the whole Ie — Ii domain. The presence of these curves is the most 
relevant difference with the weak-inhibition regime and the classic (mean-held) Wilson-Cowan model. 
Furthermore, the LP and H bifurcations that belong to the secondary branches give rise to further 
bifurcation curves (purple and light blue lines, respectively) in the Ie — 11 domain, as shown in Fig. ( |10[ ). 

Interestingly, since the branching point bifurcation increases the dimension of the network from 2 (for 
A/ < 0, see Eq. (§) to 3 (for A/ > 0, see Eq. (dll)), the network can exhibit more complex dynamics, 
such as quasi-periodic motions originated from the torus bifurcations. The biological importance of 


quasi-periodic oscillations in neural communication was discussed in 57 


Now we want to study the local bifurcations from the analytical point of view. We start by considering 
the BP bifurcations, which are dehned by the condition A/ = 0, as we saw before. From Eq. ([^ this 
condition implies: 


{ill (BP)) 


N -I 
Tl I Jll I 


( 20 ) 


SO the solutions of this equation are: 
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Figure 9: Codimension one bifurcation diagrams for fiE (left) and fj.i (right) as a function of Ie, for Jn = —100 
and 1 1 = —10. The colored curves describe the same bifurcations as in Fig. Q. Besides, we observe an LPC 
bifurcation (dark green line) and a TR bifurcation (gray line). 




Figure 10: Codimension two bifurcation diagram on the Ie—Ii plane for Ju = —34 (left) and Jii = —100 (right). 
In addition to the bifurcations already displayed in Fig. 0- we stress the presence of new ones. The branching 
points form two curves (light green lines), that define the values oi Ie — h that bound the secondary branches 
of equilibrium points (see the violet curves in Figs. (|^ and @)- The bifurcations originated on the secondary 
branches are differentiated from those originated on the primary one. Specifically, we show Andronov-Hopf and 
saddle-node curves (purple and light blue lines, respectively). In addition, we display the torus bifurcation curve 
(gray lines). 
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Now, from the second equation of the system ([^ (we can also use Eq. since for A/ = 0 they are 

equivalent) we get: 
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while from the first equation of © we get: 


( 22 ) 
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(23) 


where m/(BP) and mb (BP) are given by Eqs. (21) and (22) respectively. Since mb (BP) depends on 
//, Eq. ( [2^ defines two explicit functions Ie = IF± (Ii), that provide the curves on which we have a 
branching point bifurcation (see the light green lines in Fig. (10) for Jjj = —34 and Jjj = —100. More 
details can be found in SubSec. (S4.2.1) of the Supplementary Materials). 

The points where the H and BP curves meet each other define the ZH bifurcation. From this definition, 
we see that they can be calculated analytically from the conditions Ao,i = itw and A/ = 0, from which 
in turn we get: 
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As usual, if we substitute these expressions of the membrane potentials in Eq. ([^ or ( [II| ), we obtain the 
coordinates of the ZH points in the Ie — I I plane. 

As we said, on the secondary branches that are generated by the branching points, new bifurcations 
can occur (in the case Nj = 2, see for example the LP and H bifurcations in Figs. ([^ and Q, and the 
corresponding light blue and purple curves in Fig. ®), also new branching points (for Nj > 2), from 
which tertiary branches emerge, and so on. To study them, according to bifurcation theory, we need 
the Jacobian matrix of the network on the secondary (tertiary, and so on) branches, as we will explain 
more clearly in Sec. (3.4|. As usual, we focus on the case Nj = 2, therefore we can determine the local 
bifurcations on the secondary branches by means of Eq. (12). 

Now we start with the LP bifurcations. We know that they are defined by the condition that one of 
the eigenvalues of (12) is equal to zero. From it, as explained in SubSec. (S4.2.3) of the Supplementary 
Materials, we obtain that: 
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So if we invert Eq. (24) and we use the solution of Eq. (131, we obtain the expression of mb as a function 


of M/, 0 - If we replace the solutions mb and m/,i in the system (11), we get parametric equations for Ie,i 

as a function of a single parameter, which is now defined as D M/,o- These equations are an analytical 
description of the light blue curves shown in Fig. (10) (right) for J// = —100. 


Similarly, for the H bifurcations we obtain the condition: 
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so again it is possible to describe these bifurcations analytically, obtaining the same results we got 


numerically in Fig. (101 for Jjj = —34 and J// = —100 (see the purple curves in both the panels). 


However, unlike the primary branch, our theory does not allow us to calculate the range of the parameter 
0 = M /,0 on the secondary branches, since the resulting equations that define the range are analytically 
intractable. In the same way, now it is not possible to calculate explicitly the coordinates of the new BT 
bifurcations, where the LP and H curves that emanate from the secondary branches meet each other. 
Therefore analytical approximations or numerical schemes must be used to evaluate them, but this is 
beyond the purpose of the article. 

As we did for the weak-inhibition regime, we conclude this section by describing briefly the effect of 


the variation of the remaining synaptic weights. As we said at the end of SubSec. (3.21, all the variations 


that occur on the primary branch are qualitatively similar for different values of J//. On the other side, 
now we want to analyze the behavior of the BP curves and of the bifurcations on the secondary branches. 
For Jji = —100 and increasing Jee^ the most notable phenomenon is the overlapping between the oblique 
parts of the BP and the LP curves of the primary branch. The latter finally collapse on each other and 


split in two disjoint parts as in the case Ju = —10 discussed in SubSec. (3.2), while the bifurcations on 
the secondary branches do not show any interesting variation. Furthermore, when Jee ^ 0, we observe 
first of all the disappearance of the ZH bifurcations. This occurs because the H curves on both the 
primary and the secondary branches do not meet the BP curve anymore. If we further decrease Jee^ the 
two CP bifurcations on the LP curve of the secondary branches get closer and closer until they annihilate 
each other and the curve disappears. Clearly this phenomenon implies also the disappearance of the BT 
bifurcations on the secondary branches. For smaller values of Jee, the H curves on both the primary and 


secondary branches disappear, and finally also the LP curve on the primary branch (see SubSec. (3.2)) 
and the BP curve. To conclude, for large \ Jei\ or large Jje, the LP curve on the secondary branches 
disappears again through the annihilation of its CP points (as explained above), while on the other side, 
when at least one of the two synaptic weights is small, we do not observe any interesting variation of the 
bifurcations on the secondary branches. 


3.4 The case with generic Nj 

The same analysis can be performed on networks with a generic number Nj of inhibitory neurons. When 
A/, as given by Eq. ([^, goes to zero, we observe in general the formation of secondary branches from 
the primary one. This means that an inhibitory membrane potential becomes different from the others, 
so we can reinterpret the system as a network with an excitatory population with Ne neurons, and two 
inhibitory populations, one with one neuron, and the other with Nj — 1 neurons. Furthermore, when 
we change the current Ie (while keeping A/ > 0) with // fixed, at some point one of the eigenvalues 


of the Jacobian matrix of the system (10) (see the Thm. (SI) of the Supplementary Materials for their 
analytical calculation) may go to zero, generating a new branching point on the secondary branches. In 
this case, we observe the formation of tertiary branches, so now one of the previous Nj — 1 identical 
inhibitory membrane potentials becomes different from the others, and so on. 

In Fig. 0 we show an example of formation of secondary branches in the case Nj = 4, obtained 
as usual for the values of the parameters reported in Tab. ([^. As the reader may easily see, there is 
now an important difference compared to the case Nj = 2. For a network with only two inhibitory 


neurons, Eq. (11) implies that they both have the same codimension one bifurcation diagram (see the 
right panels of Figs. (|^ and (|^). This is just a special case, because in general, for Nj > 2, we observe a 
symmetry-breaking not only at the level of the inhibitory membrane potentials, but also at the level of the 


codimension one bifurcation diagram. Indeed, from Fig. (11) we see that one inhibitory neuron, which is 


chosen randomly by the system, has a different diagram compared to the others. This is another example 
of symmetry-breaking that occurs in the system, as an implicit consequence of Eq. (10) for Nj > 2. 


It is important to underline that, even if the BP curve (23) is defined for every Nj, numerically we 


observe the formation of new branches only for Nj even. In principle this may be proved rigorously by the 
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Figure 11: Codimension one bifurcation diagrams for Ni = 4. The left panel shows the diagram of the excitatory 
neurons, the central panel that of three of the inhibitory neurons, while the right panel shows the diagram of the 
remaining neuron in the inhibitory population. Interestingly, now the inhibitory neurons have not only different 
membrane potentials, but also different codimension one bifurcation diagrams. The neuron with the different 
diagram is chosen randomly by the system, so this is another example of symmetry-breaking that occurs in the 
network for strong inhibition. 


SO called Lyapunov-Schmidt bifurcation equation 58 , but since the proof is rather technical and beyond 


the purpose of this article, we do not report it here. 

For a generic Nj, local bifurcations can still be calculated analytically as we showed before. Indeed, 
from the second equation of (10), it is possible to express iV/ — 1 inhibitory membrane potentials as 


functions of the remaining one, which can be used as a parameter for the parametric equations in the 
codimension two bifurcation diagram. 

To conclude, we observe that now new kinds of bifurcations can appear, which do not occur in 
the case iV/ = 2. For example, for TV/ = 4, if the network has four different inhibitory potentials, 
the characteristic equation of the Jacobian matrix has the form p (A) = (A —^p//(A), where 
Pr (A) (the characteristic polynomial of the reduced Jacobian matrix introduced in Thm. (SI) of the 
Supplementary Materials) is a fifth order polynomial. This means that in principle, for some values 
of the parameters, the network may have two pairs of purely imaginary eigenvalues. This condition 
corresponds to the formation of a so called double-Hopf bifurcation, which in turn may imply a local 


birth of chaos (e.g. 59 60 ). Indeed, for A/ > 0 the dimension of the system is larger than 2 due to the 
BP bifurcations, therefore according to the Poincare-Bendixson theorem the network may exhibit chaotic 
behavior. 


3.5 Differences between our approach and the mean-field theory 

In this section we explain in detail why the mean-field theory cannot describe the branching point bifur¬ 
cation. Given the system ([^, the mean-field theory makes the assumption that within each population 
the membrane potentials are independent and identically distributed. Therefore, by hypothesis, it forbids 
the presence of heterogeneous solutions, like those that emerge from the branching point bifurcation. Due 
to this assumption, Eq. (§ can be reduced to a system of two differential equations, according to the 
mean-field theory developed by Sznitman, Tanaka, McKean and others 61-68 : 


= ':^Ve (t) + ReJee^ [-s^e (Ve (t))] + RiJei^ {Vi (t))] 4 - Ie 

= ^Vi (t) + ReJieE We {Ve (t))] + RiJnE [M {Vi (t))] -f h 


(26) 
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where Ra = lira (namely the ratio between the population a and the whole network in the thermo- 

dynamic limit, so in our case Re = 0.8 and Rj = 0.2), while Va represents any membrane potential in the 
population a. Moreover, E [•] denotes the average over trials at a given time instant, and it means that 
the system is generally supposed to be stochastic. Stochasticity can be introduced in different ways, for 
example through Brownian motions, random initial conditions, or random synaptic weights [17| . In this 
article we are considering a deterministic system, therefore we get simply E [si/ct {Va (f))] = ^ (14, (f)). In 
this way, the neural network is described by a system of two coupled equations in the unknowns Vej (t), 
whose Jacobian matrix is: 




— + ReJeejV^ (mb) RiJei^i ifj-i) 


ReJie^e (mb) 


~^ + RiJiisVl (/i/) 


(27) 


Its characteristic equation is: 


mf / \ mf \ ^ I T mf \ mf , mf r\ 

^ ^0,1 + c = 0 


where: 


a”' =1 

=- 1 - ReJee-s^e (p-e) — RiJiisVi (hi) 

TE Tl 

=-( — RiJiijVj {ill) H- ReJeesVe (mb) ) + ReRi {JeeJii — JieJei) sV'e {pe) sVj {fii) 

TeTI \Te Tl ) 


From Eq. ([7|) it easy to see that lim Aq i = A™J. The only difference between Aop and A“J is in the 

ratios that multiply the synaptic weights (or for the finite-size network, and Ra in the mean- 
field case). This difference, due to the absence of self connections, is small for large networks. So in a 
sense, when compared to a finite-size network, the mean-field approximation takes into account only the 
information provided by Ao,i, and neglects that of Xe,i- Clearly Xe is always negative, therefore it never 
affects the changes of dynamics of the system. However, in a finite-size network A/ can change sign, 
generating a branching point bifurcation. The mean-field approximation neglects this information, and 
it can be seen as a consequence of the fact that lim A/ = —1^. In other words, in the thermodynamic 

A/—foo 


limit the eigenvalue A/ is always negative, therefore it cannot generate branching point bifurcations. In 
more mathematically rigorous terms, we get that the center manifold 26 of the network is not affected 
anymore by A/ for N —>■ oo, so that the dynamics is governed only by Aop. This clearly proves that the 
mean-field approximation oversimplifies the description of the network, since it is able to describe only 
the primary branch. 


3.6 Finite-size effects are stronger for biologically plausible anatomical con¬ 
nections 

By comparing the right panels of Figs. ^ and it is clear that the magnitude of the branching point 
finite-size effects is proportional to the magnitude of J//. Indeed, from these figures we can see that. 
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for a given Ie, the difference between the membrane potentials of the primary and secondary branches 
increases with Moreover, from Fig. (lOl we observe that for J// = —100 the codimension two 


diagram is more complex than the case Ju = —34, due to the formation of LP, CP and BT bifurcations 


on the secondary branches. On the other side, in SubSec. (3.51 we have seen that the branching point 


bifurcations disappear in the thermodynamic limit. This is a consequence of the increasing number of 
incoming connections M/ = TV — 1, which implies that lim A/ = — — <0. Therefore from all these 

Af—>oo 

observations we conclude that the magnitude of the finite size effects is related to the ratio in the 
formula of A/ (see Eq. ([^). In other terms, both inhibition and the topology of the anatomical connections 
determine the strength of the finite-size effects. In particular, here we want to discuss the importance of 
the parameter M/. 

For this reason, we extend our analysis to the case of a more realistic connectivity matrix (to simplify 
matters, we consider a purely inhibitory neural network, since it is sufficient to show branching point 
bifurcations). For example, we can consider the block-circulant topology with 

circulant-band blocks introduced in 


17 


J=Jii 




(i) _ 


25(0) Q3PI 

23(1) 231^1 


25 ( 0 ) 


l-<5io 1 1 0 
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( 28 ) 


0 1 


1 0 


■ • 1 — 5i0 1 

0 1 • • • 1 1 - (5io 


where ..., are G x G circulant matrices (so that FG = N), with bandwidth -|- 1, for 

j = 0,..., E — 1. The network equipped with these synaptic connections can be interpreted as a collection 
of F neural masses with G neurons each. If we define: 


H{x) = 


X < 0 


X > 0 


then Mo = 2^0 - -ff (Co - [f J + ("I)'') is the number of connections that every neuron in a given 

mass receives from the neurons in the same mass. Furthermore, Mi '=^ 2^^ -|-1 — iJ -|- (—1)*^^, 

for i = 1,...,E — 1, is the number of connections that every neuron in the jth mass receives from 
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the neurons in the (i+j)th mod F mass, for j = 0,...,F — 1. Therefore we conclude that Mj = 


F- 1 


E F —1 
z =0 


LfJ +(-if) 


Now, if we suppose that the membrane potentials are homogeneous, the eigenvalues of the correspond¬ 
ing Jacobian matrix are: 


AmG+n — ^ 


F-1 

+ ^ F-l + Y.9iFiuG) 


Ml 


1 I Jj 


-i + 


i=Q 

F-1 


{fj,i 




i=Q 


g{n,^i,G) = < 


^( 2 Si + l) 




- 1 , 


m = 0, Vn 


s/j , m 7 ^ 0, Vn 


(29) 


'2^i-H (Ci - Lf J + (-1)'') , n = 0, VCi 
-1, n / 0, = l_f J 


n / 0, < Lf J 


with TO = 0,_F — 1 and n = 0,G — 1. They depend on the ratio where M/ does not necessarily 
diverge in the thermodynamic limit (consider for example the case when G oo for F fixed, and the 
parameters are finite and independent from G). Therefore, for N large enough, this topology exhibits 
stronger finite-size effects than the fully connected network. 

To conclude, we also underline that, according to 17 , if Mj does not diverge for iV —>• oo, the 


neurons do not become independent, therefore Sznitman’s mean-field theory cannot be used to simplify 
the description of the network. Moreover, from (29) we see that many of the eigenvalues XmG+n are 


distinct, since the reduced number of connections breaks the degeneracy of the system (compared to the 
fully-connected network, where A/ has algebraic multiplicity Nj — 1). For this reason we argue that they 
generate a multitude of branching point bifurcations, not just two as in the fully-connected case, making 
plausible connections even more interesting from the biological and mathematical point of view. 


4 Discussion 


We proved the emergence of complex dynamics in small neural circuits, characterized by strong finite-size 
effects, that cannot be accounted for by the mean-field approximation. We showed, through a detailed 
numerical and analytical analysis of the bifurcations, that small symmetric neural networks undergo 
branching point bifurcations through spontaneous symmetry-breaking, that leads to the formation of 
strongly heterogeneous membrane potentials in the inhibitory population. This result is obtained when 
we increase the strength of the synaptic weights in the inhibitory population, and clearly falsifies the 
mean-field hypothesis of identically distributed neurons. This is a very interesting feature of our model, 
since from the simple assumption of identical neurons, it is able to exhibit heterogeneous membrane 
potentials, as in real networks. 

From a biological point of view, strong inhibition may correspond to anesthetized neurons, since 
it was proved that some kinds of anesthetics, such as propofol, thiopental and isoflurane, act on 7 - 
Aminobutyric acid (GABA), the primary inhibitory neurotransmitter in the brain 69 . Our results prove 


that the dynamics repertoire of neural populations can be completely different when the brain is under 
the influence of drugs, with important consequences on experiments with anesthetized animals. They also 
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underline the importance of the synaptic weights in determining the dynamical behavior of the neural 
network. As in the theory of hallucinations of Ermentrout and Cowan 28 , the spontaneous symmetry¬ 


breaking is caused by drugs that increase the strength of the synaptic weights and generate alternative 
neural patterns. However, it is important to observe that while in their theory the network undergoes a 
symmetry-breaking through an increase of the excitatory weights, in our model symmetry is broken by 
an increase of the inhibitory ones. 

The analysis we performed can be used to understand not only the dynamics of neural masses at 
the mesoscopic scale in big animals, but also the behavior of tiny brains such as those of rotifers and 
nematodes. However, it is important to observe that in this article we restricted the bifurcation analysis 
to small neural circuits, because of our hypothesis of all-to-all synaptic connectivity between neurons. 
This hypothesis allowed us to reduce the complexity of the analytical formulas, but predicted that the 
magnitude of the finite size-effects is inversely proportional to the number of incoming connections per 
neuron, that for a fully connected topology is proportional to the network’s size. Therefore, after relaxing 
the hypothesis of all-to-all connectivity, in order to study the behavior of biologically more plausible 
networks, we argued the formation of multiple branching point bifurcations and a stronger heterogeneity 
of the membrane potentials, as a consequence of the reduced number of synaptic connections. This means 
that in principle strong finite-size effects can occur also in large networks, if their anatomical connections 
are sufficiently sparse. 

We also observe that our study may be extended to gain some insights into the oscillations that emerge 
from Hopf bifurcations. In the codimension one bifurcation diagram, when the current /_e is close to the 
one that generates a Hopf bifurcation, the amplitude of the oscillation that appears from the H point is 
small. For this reason, the solution of the system ([^ can be approximated by a truncated Fourier series. 
Semi-analytical expressions of the amplitudes of the harmonics can be obtained by applying the so called 
harmonic balance methods, which were developed especially in the context of electrical engineering 70 . 
In principle, these tools may be used to study the formation of torus bifurcations and eventually the 
transition to chaos 
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Furthermore, we stress the fact that the analytical study of local bifurcations introduced in our work 
can be easily extended to neural networks composed of several populations, by means of the Thm. (SI) 
reported in the Supplementary Materials. In this way we can shed new light on the behavior of more 
complex and biologically plausible networks, also when numerical simulations become prohibitive. So for 
example we may think to extend our model to describe the interaction between six neural populations. 
These can be thought as a coarse description of the six layers of the cerebral cortex, therefore this system 
can be interpreted a simplified model of a cortical column. In general we obtain that in a network 
with tp populations, a codimension bifurcation is described by a (fp — ‘^)-dimensional manifold, whose 
parametric equations contain tp —^ independent parameters. The only exception, as in the case with two 
populations considered in this article, is the BP manifold, whose equation can be written by expressing 
any current as an explicit function of all the others (see Eqs. (211 -|- (22) -|- (23)). So for example, in a 
network with three populations A, B, C, in the codimension three diagram spanned by the currents I a, 
Ib, Ic we get that: 


• the bifurcations LP, H etc are described by surfaces with two parameters; 

• BT, CP etc by lines with one parameter; 

• codimension three bifurcations by points; 

• the BP bifurcations are described by surfaces with explicit formulas 1a = A- {Ib,Ic)- 

However, it is important to observe that a complete classification of bifurcations with > 3 is still 
missing in the literature, due to their complexity. Nevertheless the method introduced in this article in 
principle can be used to describe analytically some of them for any 
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We underline that the Thm. (SI) can be used to describe also networks of interconnected neural masses 
(where each one is described by the system in Fig.[^. Interestingly, in the special case when the masses 
are identical, the symmetry group is more complex than that of a single mass. As we proved, we may 
observe spontaneous symmetry-breaking within a mass, through branching point bifurcations. However, 
since a network of masses has a larger symmetry, now symmetry-breaking can occur also between masses, 
which means that masses can behave differently even if they are identical. We found that in each mass 
the inhibitory membrane potentials can oscillate in synchrony or out of synchrony (on the primary and 
secondary branches, respectively). So in principle we may observe the formation of particular patterns 
of activity known as chimera states 75 , where inhibitory neurons in one mass are synchronized while 


those in the other are not, even if the two masses are described by identical neural equations. 

To conclude, we observe that our model can be extended and generalized in many other ways, through 
the introduction of delays, stochasticity, synaptic plasticity and more realistic topologies. All these add¬ 
ons allow us to improve the biological plausibility of the model, without losing the possibility to investigate 
analytically the complexity of its dynamics. 
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